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Abstract 

We introduce a general implementation of the recently proposed homogenization theory 
[Tsukerman, J. Opt. Soc. Am. B 28, 577 (2011)] allowing one to retrieve all 36 linear con- 
stitutive parameters of any 3D metamaterial with parallelepipedal unit cells. The effective 
parameters are defined directly as linear relations between pairs of coarse-grained fields, in 
contrast with methods where these parameters are obtained from reflection/transmission 
data or other indirect considerations. The method is applied to plasmonic metamaterials 
with spherical gold particles and split-ring resonators (SRR), respectively. In both cases, 
the expected physical behavior is reproduced almost perfectly, with no unphysical artifacts. 
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I. INTRODUCTION 



In recent years, very extensive research effort has been devoted to investigating 
and designing new artificial materials, called metamaterials, with staggering new 
optical properties. Some spectacular phenomena realized with metamaterials are ar- 
tificial magnetism negative refraction {2, 3], super resolution 4}], electromagnetic 
cloaking j^] and electromagnetically induced transparency {(3- 9]. Although these 
properties are not found in nature, metamaterials are in fact made of conventional 
materials that are structured on a scale smaller than the free-space wavelength Ao 
of light: the unit cell a of metamaterials is typically in the range a ~ 0.1 — 0.4Ao- 
Consequently, one may look for effective material parameters describing the optical 
properties of such nanostructured composite materials on the macroscale. The def- 
inition and computation of such effective constitutive parameters is, however, not 
a trivial task and in general involves numerical simulation, as accurate analytical 
characterization of the electromagnetic fields and parameters can only be obtained 
in special cases, e.g. for a lattice of subwavelength spheroidal particles lo|-12] or 
simplified circuit-like models of split-ring resonators (SRR) [13J]. 

The literature on homogenization of metamaterials is now so vast that we shall not 
attempt to review it here. Instead, the key distinguishing features of our approach 



are highlighted be 



available; see also |19| and [20 1. 



ow. Several monographs 14Ml6| and review papers 13, ll8| are 



In the following we would like to comment on conceptual differences between pa- 
rameter retrieval methods and first-principles derivation. The linear relation between 
the four coarse-grained electromagnetic fields E, H, D, B defined by proper averag- 
ing 2^| of the "microscopic" fields e, d, b can, in general, be represented by a 6 x 6 
matrix that is usually (and especially for plasmonic materials) frequency-dependent. 
There are two general routes for finding the entries of this matrix. In the first one 



the effective parameters of a material slab are inferred from the S'-parameters, i.e. 



from reflection and transmission of a plane wave through that slab 
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The second route of analysis involves the fundamental definition of effective 



rameters as relationships between the coarse-grained ("macroscopic") fields 
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281 ] . To this end, the coarse-grained fields must be clearly and unambiguously defined 
and linear relationships between them found. 

It is important to appreciate that the first route - parameter "retrieval" - is, 
despite its great practical significance, fundamentally different from the second one 
and should not be viewed from the same perspective. A simple analogy may help 
to clarify the difference. Molecular properties of a given chemical substance can be 
either inferred from spectroscopic measurements or, alternatively, obtained from an 
ab initio quantum-mechanical analysis. Clearly, these two approaches differ substan- 
tially - they can and should go hand-in-hand but each of them is important in its 
own right. 



The methodology proposed in 20] and implemented in this paper is akin to the 
first-principles quantum-mechanical analysis in the analogy above. These first prin- 
ciples are Maxwell's equations as well as the very definition of effective parameters as 
relations between (pairs of) coarse-grained fields inside the material. In the absence 
of this fundamental definition, it would be truly surprising that the effective parame- 
ters measured for slabs of different thickness (let alone different shapes) would come 
out in agreement, even for natural materials in classical electrodynamics. The very 
consistency of effective parameters derived from different measurements and different 
geometries stems from the underlying fundamental physical characteristic - namely, 
again, effective parameters as relations between coarse-grained fields. 

Thus one distinguishing feature of our approach is that it is direct - the effective 
parameters are not inferred from some ancillary data but are derived from the fun- 
damental definition. This way, all 36 material parameters are properly and uniquely 



defined and can be found. In contrast, the "retrieval" procedures derive the effective 
parameters indirectly and allow one to determine only a subset of these parameters. 
Being in essence an inverse problem, parameter retrieval may be ill-posed and may 
suffer from the m ultip licity of solutions due to the branch ambiguities that must be 



properly handled [26]. 



Another key distinguishing feature of our methodology has to do with the way the 
coarse-grained fields are defined 2(3] . In a significant, but well justified, departure 
from the traditional ways of field averaging, we require that the coarse-grained fields 
be constructed to satisfy exactly Maxwell's boundary conditions at all interfaces - 
namely, the tangential continuity of the E and H fields and the normal continuity 
of D and B fields. Violation of these boundary conditions would be equivalent 
to spurious surface currents and/or charges, with the related nonphysical artifacts. 
For metamaterials, standard volume averaging (or, more generally, averaging by 
convolution with a smooth Gaussian-like mollifier) is not adequate in this regard. 
Indeed, if such averaging is applied within each of two abutting materials separately, 
then nonphysical field jumps will generally occur at the interface. If, alternatively, 
the moving average is applied across the interface, the fields will be smeared in a 
boundary layer, which is again nonphysical. Such spurious effects are negligible when 
the lattice cell size is vanishingly small relative to the wavelength (which is usually 
the case for natural materials), but for metamaterials nonphysical conclusions may 
result because the cell sizes are appreciable. 

Thus the microscopic fields, in our approach, are subject to two different types of 
interpolation. The E and H fields are obtained from e and b/^ by an interpolation 
that preserves tangential continuity, while D and B must be obtained from d and 
b by an interpolation that preserves normal continuity [20.]. The fact that h (i.e., 
b//xo) and b must be subject to different types of interpolation explains, among 
other things, why the macroscopic H and B fields may be different even though 



their microscopic counterparts are the same (for intrinsical 
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y nonmagnetic media); 



30). 



this paradox is at the heart of "artificial magnetism" 

Although the theory of [2GJ, reviewed in the following section, uniquely defines 
the full 6x6 parameter matrix, it has so far been applied to several test cases 
involving a subset of the 36 parameters [20]. This paper describes a general im- 
plementation of this method in 3D for a parallelepipedal unit cell, for any type of 
inclusions and without any a priori assumptions about the effective parameters. The 
implementation is described in Section [II] Sec. IHIl is devoted to instructive test cases 
of plasmonic metamaterials with gold spherical particles and split-ring resonators 
(SRR). The metamaterial with a lattice of spheres was chosen because Lewin's the- 
ory 12] is available for comparison. The SRR-metamaterial is interesting because 



of its ubiquity and the bianisotropic response [22|, |3jJ. In Sec. [TV] the computed 
parameters of the SRR-metamaterial are compared with the parameters obtained by 
the S-parameter methods. Conclusions are given in Sec. [V] 



II. THE METHOD AND ITS IMPLEMENTATION 

A. Vectorial Interpolation for the Coarse-Grained Fields 

In this section, we implement the direct definition of parameters outlined in the 
Introduction. Namely, the coarse-grained fields are rigorously defined and a linear 
relationship between them established. As already emphasized, the interpolation 
procedures for constructing the "macroscopic" fields must be chosen carefully, so 
that the normal or tangential continuity of the respective fields across all interfaces 
is honored. 

Even though the relevant theory was already presented in it is revisited 
here, to make the paper more easily readable and self-contained. Some details and 
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Figure 1. a) Sketch of the unit cell of size 2a x 2b x 2c containing an arbitrary metallic 
inclusion, b) Cross section of a split-ring resonator (SRR) lying in the xy-plane with a 
bending radius and gap of 35 nm and 9 nm, respectively. The width and thickness of the 
SRR are also 35 nm. 

implementation issues that have not been previously discussed are also presented 
here. We assume that the unit cells are parallelepipeds containing arbitrary inclusions 
[Fig. Ufa)] with a linear complex dielectric permittivity e(r,oj). The size of the unit 
cell is 2a x 2b x 2c along the x, y, z directions, respectively. 

The interpolation procedure that preserves tangential continuity is effected by 
vectorial interpolation functions like the one shown in Fig. [21 in a 2D rendition for 
simplicity. The circulation of this function is equal to one along one edge (in the 
figure, the vertical edge shared by two adjacent lattice cells) and zero along all other 
edges of the lattice. Within each cell, the interpolation function varies linearly. A 



formal expression for four of these functions is 



wi_4 = ^ ( i ± 7) fi ± - ) x. (i: 



k \ bJ \ Cj 

Another eight functions of this kind are obtained by the simultaneous cyclic permu- 
tation of (x, y, z) and (a, b, c) in the expression above. For each lattice cell, there are 
12 such interpolating functions altogether (one per edge). Each function w a has unit 
circulation along edge a (a = 1, 2, . . . 12) and zero circulations along all other edges. 
The coarse-grained E and H fields can then be represented by interpolation from 

the edges into the volume of the cell as follows: 

12 12 
E = J^[e] a w Q , /x H = ^[b] a w Q , (2) 

where [e] a = J e-dl is the circulation of the (microscopic) e field along edge a; similar 
for the b field circulation [b] Q . In the calculation of the circulations, integration along 
the edge is always assumed to be in the positive direction of the respective coordinate 
axes. 

We now move on to the second kind of interpolation that preserves the normal 
continuity and produces the D, B fields from d and b. A typical interpolating func- 
tion (2D rendition again for simplicity) is shown in Fig. [3j The flux of this function 
through a face shared by two adjacent cells is equal to one; the flux through all other 
faces is zero. There are six such functions per cell: These vector functions can be 
formally expressed as 

v " = {ss( 1:t i)*'^( 1:t ?)''8b( 1:t f)*}' (3) 

and can be used to define the coarse-grained D and B fields by interpolation from 
the six faces into the volume of the unit cell: 

D = ^[[d]]^, B = f^MfiVp, ( 4 ) 



Figure 2. (From [2JJ].) A 2D analog of the vectorial interpolation function w a (in this 
case, associated with the central vertical edge shared by two adjacent cells). Tangential 
continuity of this function is evident from the arrow plot; its circulation is equal to one 
over the central edge and to zero over all other edges. 



where [[d]]^ = f„ d • dS is the flux of d through face (3 ((3 = 1, 2, . . . , 6); similar for 
the b field. In the calculation of fluxes, it is convenient to take the normal to any 
face in the positive direction of the respective coordinate axis (rather than in the 
outward direction). 

The coarse-grained E and H fields so defined have 12 degrees of freedom in any 
given lattice cell. From the mathematical perspective, these fields lie in the 12- 



dimensional functional space spanned by functions w„; we shal 



keep the notation 



W CUT \ of [20| for this space: the 'W honors Whitney [32| (see (20| for details and 
further references) and 'curl' indicates fields whose curl is a regular function rather 
than a general distribution. This implies, in physical terms, the absence of equivalent 
surface currents and the tangential continuity of the fields involved. 

Similarly, D and B within any lattice cell lie in the six- dimensional functional 
space Wdi v spanned by functions v^. Importantly, it can be shown that the div- and 
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Figure 3. (From [20].) A 2D analog of the vectorial interpolation function vp (in this case, 
associated with the central vertical edge). Normal continuity of this function is evident 
from the arrow plot; its flux is equal to one over the central edge and zero over all other 
edges. 

curl-spaces are compatible in the following sense: 



That is, the curl of any function from W CUI \ (i.e. the curl of any coarse-grained field 
E or H denned by (j2J)) lies in W& v . Because of this compatibility of interpolations, 
the coarse-grained fields, as proved in [20|, satisfy Maxwell's equations exactly. By 
construction, they also satisfy the proper continuity conditions at all interfaces. Thus 
the coarse-grained fields are completely physical. 

B. Approximation of the Fields and Linear Relationships 

The vectorial interpolations above are written for the exact fields that are, obvi- 
ously, in general unknown. One therefore looks for a good approximation via a linear 
combination of basis modes ip^ 20]: 



V x W cmi G Wdiv 



(5) 
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\I/ and all ^ 7 denote, in the most general case, six- component vectors comprising 
both microscopic fields; e.g. ^ eh = {^ e ,^ h }, etc.; complex coefficients. 

To each basis wave tp y , there correspond the curl-interpolants E 7 (r) = W CU ri([V ; 7]i-i2 
and H 7 (r) = VV cur i([-0 7 ] i_i2). To define these interpolants (coarse-grained fields), 
one computes the edge circulations of each basis mode and applies the interpolation 
formula ([2]). Similarly, there are the div-interpolants D 7 (r) = Wdi v ([[^ 7 ]]i-6) and 
B 7 (r) = Wdiv([[V ; 7 ]]i-6)) obtained by computing the face fluxes of the basis modes 
and then applying interpolation (0J. 

For all basis waves a at any given point r in space, we seek a linear relation 

< B (r) = 0(r)Vf H (r). 

The 6x6 parameter matrix Cl characterizes, in general, anisotropic material behavior 
with (if the off-diagonal blocks are nonzero) magnetoelectric coupling. Similar to i/j^, 
the six-component vectors ip^ B comprise both fields, but coarse-grained; same for 
. In matrix form, the above equations are 

^ DB (v) = n(r)^ EH (r) (6) 

where each column of the matrices ty DB and l M contains the coarse-grained fields 
of the respective basis function. 

In this paper, the basis modes ^ 7 (7 = 1,2,..., 12) are chosen as a set of twelve 
Bloch waves propagating in six coordinate directions (±x, ±y, ±2), with two different 
polarizations per direction. This leads to an overdetermined system of 12 constitu- 
tive relations for the 6x6 parameter matrix. This system is solved in the least- 
squares sense. Qualitatively, a good least-squares fit implies that different modes 
in the metamaterial can be described accurately enough via the same set of effec- 
tive parameters, i.e. the material is "homogenizable" . A poor fit corresponds to a 
"non-homogenizable" material; this can also be viewed as strong "spatial dispersion" . 



(Rigorous quantitative definitions to be given in [33| ) . 
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Mathematically, the material parameter matrix is appropriately expressed via the 



pseudoinverse 



34): 



Q(r) = ^ DB (r) (y EH (r)) . (7) 

This expression is equivalent to solving the least-squares problem. 

The parameter matrix Cl(r), being a function of position r, is an interim quan- 
tity. The effective parameters are found by averaging 0(r) over the whole unit cell 







\B(p,«)j 






(8) 



201 ] . This accomplishes our goal of finding an effective 6x6 material matrix that 
describes, in the frequency domain, linear constitutive relations between the different 
electromagnetic fields: 

'Hp) 

The frequency-dependent parameter matrix Q comprises four 3x3 matrices with e 
and p, describing the electric and magnetic response, respectively, whereas £ and £ 
characterize the magnetoelectric coupling. 

When accurate analytical approximations of the microscopic fields are available, 
they may be used to generate the basis set; however, this is feasible only in the 
simplest cases. Generally, numerical techniques have to be employed. In this paper, 
the finite element method (FEM) is used to find the basis modes in the cell. As noted 
above, these modes are chosen as twelve Bloch waves of the form e = e per exp^fc^n^r- 
r) (with the phasor convention exp[— iut}). Here e per is a periodic amplitude function, 
k B is the complex Bloch wavenumber and n<jir is the unit vector of propagation 
direction. If the Bloch wave solution is substituted into the vector wave equation, 
one obtains an eigenvalue equation for e per 

(V + ik B n dir ) x (V + ik B n dir ) x e per - k 2 ee per = 0, (9) 

where k B is an eigenvalue and e per is an eigenmode. In this paper, this eigenvalue 
solver has been implemented in the commercial FEM software Comsol Multiphysics 
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351 ] . Since the eigenproblem is written for the periodic factor e per rather than the 



total field, periodic boundary conditions are applied on all three pairs of faces of the 
unit cell [36] . By considering six directions of propagation n^r and two polarizations 
per direction, one constructs a 12-mode basis for the homogenization procedure. 
The construction of the 12-mode basis requires for each frequency three independent 
eigenvalue calculations corresponding to n<&r being chosen along the three coordinate 
axes. The reason why only three calculations are sufficient (and not six) is because 
the eigenvalue solver computes the ±kg solutions, i.e., each choice of n^r covers the 
two propagation directions ±rw. As an example, if hdir = (1,0,0) the basis modes 
for ±x-propagation with y- and z-polarization are obtained (four basis modes are 
determined per h dir ). Note that, in addition to propagating modes, the calculation 
may reveal some evanescent ones. These modes do not affect the optical properties of 
voluminous bodies but may be considered in future studies where surface waves are 
of interest. It goes without saying that, when the edge circulations and face fluxes 
are calculated for the construction of the coarse grained fields, the microscopic field 
is not, e.g., e per but rather e per exp^sn^ • r). 



III. TEST CASES 

This section is divided into three subsections. The first one discusses the seemingly 
trivial, and yet instructive, case of a homogeneous cell. In the second subsection, all 
36 effective parameters are computed for a metamaterial with spherical gold nanopar- 
ticles; to test the methodology, no a priori symmetry assumptions are built into the 
procedure. In the third subsection, a metamaterial with gold SRRs is considered and 
all 36 parameters are again found without any preconditions. In all our simulations, 
the Johnson & Christy data for the dielectric function of gold are used. 
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A. Case 1: The homogeneous cell revisited 



Effective parameters of a homogenous cell, as a "sanity check" of the proposed 
method, were already discussed in [20]. It was shown that the method indeed yields 
the exact results, i.e., e and p, are diagonal matrices with the values equal to the 
intrinsic permittivity and permeability of the homogeneous material (e and /x, re- 
spectively), and £ and ( are zero matrices (no magnetoelectric coupling). Although 
this result appears trivial, we note that in some existing theories spurious Bloch-like 
factors do appear in the effective parameters and then need to be discarded by fiat. 

To elaborate, the homogenization procedure outlined above yields for a uniform 
cell, after a direct analytical calculation, the effective permittivity 

arctan(tan(A;a)) 
en = e — > k = io^/Jie (10) 

which results in en = e if ka < 7r/2; the subscript ii indicates that only diagonal 
elements are different from zero. Thus the method is applicable when the wavelength 
is greater than double the unit cell size (4a < A). This is perfectly consistent with 
the fact that homogenization makes sense only when the unit cell is substantially 
subwavelength. 



B. Case 2: Spherical gold particle 

Effective parameters of a material with subwavelength-sized spherical particles 
are described well by the Maxwell-Garnett formula 10J, [ll| and the Lewin theory 



for a cubic lattice of spheres 



12]. This is therefore an excellent test case, where the 



new procedure can be tested against the known analytical results and the dispersion 
around the plasmon resonance wavelength is very strong, a property that any useful 
homogenization method must capture. 
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In the following, a spherical gold particle with radius tq = 20 nm in a cubic unit 
cell of half-size a = b = c = 40 nm is considered. The host material is assumed 
to be non-polarizable (or free space). As the volume fraction of the spheres is only 
Vf ~ 0.07 and the wavelength interval of consideration is 300 — 900 nm, the calculated 
results can appropriately be compared with the Lewin expression for the diagonal 
entries of the effective permittivity matrix I2I ]: 

e « = €h ( 1 + m £k ~ ) ' 

V F(6)-e h /e s U fJ 

where F{9) = 2 (sin 6 — 6 cos#)/[(6> 2 — 1) sin# + 6 cos 9], 8 = k r ^e s fi s / (e /io), e s and 
Us the permittivity and permeability of the spheroidal inclusions, respectively, and 
eh the permittivity of the host material. For gold spheres in a vacuum, e s = 6^(0;), 
fjL s = A*q and = eo- Subscript ii in Eq. (TTTj) refers to either xx, yy or zz. The 
permittivity tensor is obviously diagonal; however, this is not explicitly built into the 
homogenization procedure as an a priori assumption but rather is expected to emerge 
naturally as a result. Also, the Lewin theory does not include any magnetoelectric 
coupling (£ and ( are zero matrices), whereas the magnetic response is described as 
in Eq. (fTTj) with the substitution of e with /i. 

The procedure of [20] and Sec. [IT] lias yielded all 36 parameters that agree quan- 
titatively well with the Lewin theory (Fig. H]). In this figure, only the diagonal 
elements of the permittivity and permeability matrices are displayed. All other com- 
puted parameters are negligibly small, although not identically zero due to numerical 
errors. Specifically, the off-diagonal elements of e and p, are 10 4 — 10 6 times smaller 
than the diagonal elements, whereas the dimensionless values of c £ and c ( (c is 
the speed of light in a vacuum) are smaller by a factor of ~ 10 5 compared to the 
values of e/e and /2///o, thereby confirming that the effect of the magnetoelectric 
coupling is negligible. 
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Figure 4. The diagonal elements of a) the effective permittivity matrix and b) permeability 
matrix for a metamaterial consisting of cubic-ordered gold spheres with diameter of 40nm 
embedded in a vacuum. The unit cell size is (80nm) 3 . The ovals and arrows identify which 
curves belong to which axes. 



The computed diagonal elements of the permittivity matrix are all equal and also 
consistent with the Lewin theory [Fig. H(a)], exhibiting only a small discrepancy at 
shorter wavelengths. Even the plasmon resonance at ~ 520 nm, resulting in a strong 
dispersion and increased absorption, is well captured in the effective parameters. The 
higher discrepancy at shorter wavelengths is a natural consequence of approaching 
the homogenization limit (double the cell size) at approximately A ~ 160 nm. The 
permeability matrix [Fig. Hl^b)] is also consistent with the expected behavior; its 
diagonal elements are predominantly real and equal to ~ 1. 
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The agreement with the Lewin theory is, however, not perfect. At shorter wave- 
lengths, the real part of the permeability is smaller than Lewin's value of ~ 1, while 
the small imaginary part is slightly overestimated as compared to Lewin's theory. 
These discrepancies are a result of leaving the subwavelength regime. One should 
note that although the relative deviation is larger for nn compared to en, the absolute 
difference is approximately ~ 0.1 for both En and [in. It is also worth noting that the 
retrieved permeability is affected to some extent by the electric plasmon resonance 
and exhibits a weak perturbation near the resonance wavelength. We believe this is 
a consequence of numerical errors. For instance, the b field is obtained by numerical 
differentiation from the calculated electric field in the eigenvalue problem [see Eq. 

which affects the accuracy of the respective edge circulations and face fluxes. 
Importantly, though, both e and Jx satisfy the passivity requirement Im[£jj] > and 
Im[//jj] > 0. 



C. Case 3: Split-ring resonators 

The split-ring resonator (SRR) is the ubiquitous subwavelength meta-atom in 
metamaterial design tlfl, giving rise to artificial magnetism or negative refractive 

n 

index [21]. Even though SRR is generally utilized as a magnetic meta-atom, it also 
has, as noted in e.g. [22I, 0]> a magnetoelectric response characteristic of optical 
activity. In the following, a gold C-shaped SRR in an otherwise empty cubic unit 
cell with a = b = c = lOOnm is analyzed. Having a width and thickness of 35 nm, 
the SRR lies in the xy-plane parameterized by a radius of 35 nm and a gap of 9 nm 
[Fig. QJb)]. With this configuration, the magnetoelectric response arises due to an 
x-polarized incident wave inducing a magnetic dipole moment along the z-axis, and, 
conversely, a magnetic field along the z-axis generating an electric dipole moment 

16 



along the x-axis. Thus the effective parameter matrices take on the following form: 



( e xx ^ 



Eyy 
e zz 

o o £, N 






/ 



^ 1 ^ 
1 

L \l zz , 

^ 0^ 


\tz* o oy 



;i2) 



where ^ = — Czx due to reciprocity 



131 ] . As in the previous test of the new method- 



ology, no advance information about the symmetry and structure of the parameter 
matrices was explicitly utilized. However, as in the previous case, only the dominant 
parameters are shown in Fig. [5j 

The off-diagonal elements in e and p are 10 4 — 10 6 times smaller than the diagonal 
elements, while the elements in £ and ( are smaller than £ xz and ( zx by a factor of 
10 3 - 10 5 . 

It is evident from Fig. [3(a) that e xx exhibits a plasmonic resonance response at 
~ 1260 nm. This is due to the induced electric dipole moment pointing along the 
gap, in the x-direction. In contrast, e yy and e zz are almost frequency-indpendent, 
representing just a perturbation of the unit parameters of the empty cell by the inclu- 
sion of gold SRRs. Concurrently with the electric response, the SRR-metamaterial 
also displays a magnetic response in [i zz [Fig. Efb)], dramatically increasing the 
imaginary part while the real part shows diamagnetic and paramagnetic behavior 
below and above the resonance wavelength, respectively. This is all in accordance 
with analytical approximations of SRR by an equivalent circuit diagram [la ]. Again, 
it is noticeable how the non-resonant diagonal elements decrease slightly {\iaj \i§ < 1) 
for shorter wavelengths (limit of homogenization is at ~ 400 nm). 

For easier comparison of the the dimensionless magnetoelectric coupling param- 
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Figure 5. The dominant elements of a) the effective permittivity matrix, b) effective perme- 
ability matrix, and c) the magnetoelectric coupling matrices for a metamaterial consisting 
of cubic-ordered gold SRRs of Fig. [TJ The unit cell size is (200nm) 3 . The ovals and arrows 
identify which curves belong to which axes. 
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eters, Fig. [5](c) shows cq^ xz and — cq( zx . It is evident that £ xz and — ( zx are very- 
similar but not identical; there are noticeable deviations from reciprocity near the 
resonance. We believe that the small difference in the parameters, cq\^ xz + ( zx \, in- 
dicates the approximation error inherent in the procedure. At the resonance, where 
the difference is largest (c |^ X2 + ( zx \ — 0.1), we expect the approximation error to 
be higher. 



IV. COMPARISON WITH THE S-PARAMETER METHOD 



Even though the computed parameters exhibit the expected behavior for both the 
sphere- and SRR-metamaterial, there is no analytical result in the plasmonic regime 
for the SRR-metamaterial that can be used for direct comparison (only approxi- 



mations derived from electric circuit models are available [13|, [16|, |29_|). In order to 
quantitatively verify the calculated effective parameters from Sec. II II C\ this section 
studies and compares parameters of the SRR-metamaterial obtained with the well 



established S-parameter method originally proposed for isotropic metamaterials 



2l| 



and then extended to anisotropic and bianisotropic metamaterials of the SRR-type 



22] . In this paper, a simplified retrieval procedure of Ref . [23] is used. In these pro- 
cedures, regardless of their specific variants, the effective parameters are data-fitted 
to the complex reflection and transmission coefficients of a finite-sized slab of the 
metamaterial. In our calculations, the SRR-slab has the thickness of three unit cells 
(600 nm). In the following we denote the S-parameter method 'S' or the 'S-method'. 

The method studied in this paper will be denoted 'D' or 'D-method' as an abbre- 
viation of 'Direct', thereby emphasizing the fact that effective parameters are derived 
from the linear relations between coarse-grained fields. To bring the conventions in 
line with the S-parameter approach, in which only a subset of the effective parame- 
ters can be retrieved, the magnetoelectric coupling is described by a single parameter 
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Figure 6. Comparison of the parameters e xx , [i zz and xo for the SRR-metamaterial obtained 
by the direct method ('D') presented in this paper and the S-parameter method ('S') in 
Ref. [J. 



Xo in which £ xz = —ixo and ( zx = ixo- Hence, for the D-method the parameter is 
defined as Xo = i{£xz - C**)/2. 

The SRR-metamaterial is described by five effective parameters in the S-method, 
but for comparison it is only interesting to study e xx , fi zz and Xch as they are the 
only parameters affected by the plasmon resonance. The parameters are displayed 
in Fig El The electric response [Fig. |6^a)-(b)] is very similar for both methods, with 
only a small variation near the resonance. In contrast, the difference is substantial 
in the magnetic response [Fig. [6]^c)-(d)] with Im[/i 22 ] exhibiting a difference of a 
factor of ~ 2.5. Similarly, Re[/i 22 ] has noticeable discrepancies, with both methods 
showing a slight decrease in /i 22 at shorter wavelengths. The situation is somewhat 
different for the magnetoelectric coupling [Fig. Ete)-(f)], where the two methods 
show curves close to one another that from physically intuition behave as expected 
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with a response of Re[coXo] situated around Re[coXo] = 0. 

It is evident from the comparison above that both methods qualitatively agree 
well, but the specific values differ in some instances. In relation to the transmission 
and reflection coefficients, these parameter differences might conceivably add up or 
perhaps counterbalance each other, thereby leading to less or more similar spectra. 
To investigate this further, the reflection and transmission coefficients (r± and t) for a 
I/mi thick SRR-metamaterial slab are calculated for a normally incident axpolarized 



plane wave propagating in the ^-direction (Fig. [7J using the expressions [23] 



= 2i sm(nfc Qr» 2 + (x[;±^L) 2 ] , s 

r± [(tiz + n) 2 + (x r ) 2 }e- mk ° l - i(K z ~ n) 2 + (xj) 2 ]e**»' 1 ' 

t = (14) 

[{Hl z + n) 2 + ( X rf}e-^ol - _ n) 2 + ^2 ]e ink l ' 



where I is the thickness of the slab and n = \/ £ r xx H r zz — (Xo) 2 ^ ne refractive index 
of the metamaterial. The superscript r is used to emphasize that the effective pa- 
rameters are the relative dimensionless values {\f zz = /i 2z //io etc.). Similarly, the 
subscript ± on the reflection coefficient indicates that the reflection is sensitive to 
the indirection of incidence due to strong bianisotropy for the chosen direction and 
polarization. One should also note that the effective parameters of the S-method are 
in fact retrieved from the transmission and reflection coefficients in Eq. (fl3|) - fTH|) 
for a 600-nm thick SRR-metamaterial slab. For this reason, the spectra with pa- 
rameters from the S-method are expected to be in (almost) perfect agreement with 
the full-wave simulations of the transmission and reflection spectra. This is indeed 
also the evident from Fig. [71 where the curves from the S-method and 

full-wave simulation are coinciding. Consequently, it can be deduced that a 3-layer 
SRR-metamaterial is sufficient to retrieve the "bulk" effective parameters, i.e., the 
effective parameters of the S-method (Fig. [6]) are converged. 

The comparison of the D-method spectra with the full- wave simulations are equiv- 
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Figure 7. Comparison of the complex reflection and transmission coefficients from a 1/im 
thick SRR-metamaterial slab for the two different homogenization methods. The incident 
wave is x-polarized and propagating along the y-direction. r + and r_ correspond to the 
reflection coefficients for propagation along +y and -y, respectively. The legend 'Sim' 
denotes full-wave simulation of a 5-layer SRR-metamaterial. 

alent to comparing the two homogenization approaches. Comparing the reflection 
spectra obtained by the two methods [Fig. [7(a)- (d)], it is remarkable how well the 
curves are consistent with one another, even though the retrieved permeability \i zz 
is significantly different [Fig. [6(c)- (d)]. The difference in effective parameters seems 
to (partially) counterbalance each other. That said, it is still clear that the reflection 
weakly depends on the choice of effective parameters; this dependence is, however, 
much less rigid for transmission. The transmission coefficient [Fig. [7(e)- (f)] is only 
marginally affected by the choice of the method, with the two methods showing an 
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almost perfect agreement in both the amplitude and phase. It should be noted that 
such an insensitivity to the transmission coefficient is in fact a problem for the S- 
method, as small numerical errors in the S-parameters may substantially influence 
the retrieved effective parameters. 

A similar comparison of effective parameters for the SRR-metamaterial have also 
been conducted for the SRR-metamaterial with a smaller unit cell half-size a = b = 
c = 80nm (not shown). Naturally, the dispersion of the effective parameters is in this 
case higher, but the differences in the parameters between the two methods follow 
the same trend as for the larger unit cell (Fig. [6]). Similarly, the effective parameters 
of the D-method accurately describe the optical properties of the SRR-metamaterial 
showing only small deviations in the transmission and reflection spectra compared 
to the full wave simulations. 



V. CONCLUSION 



The paper describes a general three dimensional implementation of the homog- 



enization theory developed in 



20]. In contrast with the common "retrieval" pro- 
cedures where the effective medium parameters are inferred from transmission and 
reflection data, the proposed approach follows directly from the definition of ma- 
terial parameters as linear relations between the coarse-grained fields. To put it 
differently, our theory is aimed at understanding the "inner workings" of the meta- 
material, whereas S-parameters characterize it as a black box from the point of view 
of an external observer. These two categories of methods - u ab initio" analysis vs. 
indirect inference - complement one another. However, since these categories are 
philosophically different, their comparison should go beyond pure utilitarian advan- 
tages that one of them may or may not have in a particular case at a particular stage 
of development. 
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The key physical insight in our theory is that the E and H fields must satisfy 
the interface continuity conditions different from that of the D and B fields, and 
hence two different interpolation procedures are called for. These interpolations are 



effected by two types of linear vector functions (Section [II Bl and The fact that 

the H and B fields are subject to different kinds of interpolation explains the paradox 
of "artificial magnetism:" why the macroscopic B and H fields can be different even 
though the underlying microscopic field is the same. 

The method and its implementation do not require any a priori assumptions about 
the 6x6 constitutive matrix and allow one to find all 36 constitutive parameters. 
From the computational perspective, the procedure involves, at any given frequency, 
the calculation of 12 Bloch modes whose fields on the faces and edges of the unit 
cell are used to construct the coarse-grained fields by the two types of vectorial 
interpolation. The result is unambiguously defined coarse-grained fields obeying 
Maxwell's equations and the continuity conditions. 

Calculations for a homogeneous cell show that the method produces the exact 
result, with no spurious factors, as long as the cell size does not exceed one half of 
the wavelength. This condition must always be satisfied in practice, for the homog- 
enization to be meaningful. 

The homogenization of a metamaterial with spherical gold particles has been 
carried out and compared with the Lewin theory showing, in general, an almost 
perfect match. It is only at shorter wavelengths (the unit cell larger than one quarter 
of the free space wavelength) that deviations occur and the permeability is estimated 
to be slightly below unity. 

The effective parameters computed for a periodic SRR-metamaterial exhibit the 
expected physical bianisotropic behavior with magnetic and magnetoelectric re- 
sponses. Furthermore, the parameters are compared with S-parameter retrieval 
exhibiting a qualitatively similar behavior but some quantitative differences. It was 
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then demonstrated how the differences in parameters influence the reflection and 
transmission coefficients of an SRR-metamaterial slab. It turns out that the new 
homogenization method describes reflection and transmission very accurately. 

We believe that our successful implementation of the method will be useful in 
the future studies and design of complex nanostructures, especially if the specific 
properties of the effective material matrix are not known. Just like the S-parameter 
the new method and its implementation may also be transferred to 



procedure 
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the field of acoustics and applied to acoustic metamaterials. 
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